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ABSTRACT 

In cold dark matter (CDM) cosmogonies, low-mass objects play an important role 
in the evolution of the universe. Not only are they the first luminous objects to shed 
light in a previously dark universe, but, if their formation is not inhibited by their own 
feedback, they dominate the galaxy mass function until redshift z ~ 5. In this paper we 
present and discuss the implementation of a 3D cosmological code that includes most 
of the needed physics to simulate the formation and evolution of the first galaxies with 
a self-consistent treatment of radiative feedback. The simulation includes continuum 
radiative transfer using the "Optically Thin Variable Eddington Tensor" (OTVET) 
approximation and line-radiative transfer in the H2 Lyman- Werner bands of the back- 
ground UV radiation. We include detailed chemistry for II2 formation/destruction, 
molecular and atomic cooling/heating processes, ionization by secondary electrons, and 
heating by Lya resonant scattering. 

We find that the first galaxies ( "small-halos" ) are characterized by bursting star 
formation, self-regulated by a feedback process that acts on cosmological scales. The 
mass in stars produced by these objects can exceed the mass in stars produced by 
normal galaxies; therefore, their impact on cosmic evolution cannot be neglected. The 
main focus of this paper is on the methodology of the simulations, and we only briefly 
introduce some of the results. An extensive discussion of the results and the nature of 
the feedback mechanism are the focus of a companion paper. 



Subject headings: Cosmology: early universe — cosmology: theory — galaxies: dwarf — 
galaxies: evolution — galaxies: formation — galaxies: high-redshift — intergalactic medium — 
methods: numerical 
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1. Introduction 

In the last five years, precise observations of the cosmic microwave background (CMB), Lya 
forest, galaxy clusters, and large scale structure have produced a large amount of data that strongly 
constrain the standard cosmological model. Until now, despite some unsolved problems (Moore 
1994; Klypin et al. 1999; Navarro & Stcinmctz 2000), cold dark matter (CDM) cosmogonies remain 
the most predictive cosmological models. In CDM cosmogonies, the large galaxies that we observe 
today arc the result of the mergers of a large number of small and faint galaxies (dwarf galaxies). 
The number of dwarf galaxies observed in the local universe is larger than the number of normal 
galaxies but, today, they account for only a small fraction of the mass in collapsed objects. About 
12 billion years ago, galaxies with the mass of our Milky Way were extremely rare, and the universe 
was filled with objects similar to dwarf spheroidal galaxies. 

How small can a galaxy be? Or how small can a protogalaxy be in order to be able to form 
stars? A necessary condition for star formation is that the gas in the protogalaxy must cool to a 
temperature lower than the virial temperature of the halo in a Hubble time. If the dark matter (DM) 
mass is Mdm ~ 10^ M©, the gas in the protogalaxy (a so-called "large-halo" protogalaxy) cools 
by Hi emission lines because the virial temperature of the halo is T^ir ^ 10"^ K. But if M^m ~ 10^ 
M0 and the protogalaxy is of primordial composition, the gas cannot cool {T^ij- ^ 10^ K) unless 
molecular hydrogen (H2) is present. We call such an object a "small-halo" protogalaxy. During the 
virialization process, the gas is partially ionized {xe ~ 10~^) and heated to the virial temperature 
by shocks. Free electrons can be captured by neutral hydrogen (Hi ) to produce H^, which is 
the main catalyst for H2 formation in a gas of primordial composition (Lepp & Shull 1984). The 
chemical reaction Hi -|-H~ ^ H2 -|-e^ produces enough H2 {xh2 ~ 10~^) to allow the formation of 
galaxies as small as Mdm ~ 10^ Mq at z ~ 30, i.e., with virial temperature T^ir ~ 360 K (Tegmark 
et al. 1997; Abel, Bryan, & Norman 2000). 

"Small-halo" objects are thought to be the first galaxies formed in the universe and, if their 
formation is not inhibited by feedback, they should account for the bulk of mass in stars until 
reionization at 2; ~ 6. In Figure 1 we show the fraction of collapsed (virializcd) DM as a function 
of the DM halo mass at z = 24, 19, 16, 9, 5, 3, calculated with the Press-Schcchtcr formula (Press 
& Schechter 1974). The inserted table shows the fraction of the collapsed mass in "small-halo" 
objects {Tyir < lO'* K), normal galaxies ("large-halo" objects with 10^ K < Tyir < 10^ K), and 
clusters (T^^^ > 10^ K). The last column in the table shows the ratio of the collapsed mass of the 
"small-halo" to the "large-halo" objects. The hexagons, squares, and triangles show the mass of 1 
(T, 2 a, and 3 a perturbations of the initial density field, respectively. "Small-halo" objects are the 
dominant fraction of the DM collapsed mass before redshift z ^ 9, and they remain an important 
fraction of the DM collapsed mass until redshift z ~ 5. 

Competing feedback effects determine the fate of "small-halo" objects. Radiative feedback reg- 
ulates the formation and destruction of H2 in the protogalaxies. Mechanical and thermal feedback 
from supernova (SN) explosions could blow away the interstellar medium of "small-halo" objects. 
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Fig. 1. — Fraction of collapsed (virialized) DM as a function of the DM halo mass at z = 25, calcu- 
lated using the Press-Schechter formalism. The inserted table shows the fraction of the collapsed 
mass in "small-halo" protogalaxies {Tmr < 10^ K), normal galaxies ("large-halo" protogalaxies 
with lO^K < Tyir < 10^ K), and clusters {T^r > 10® K). The last column in the table shows the 
ratio of the collapsed mass of "small-halo" to "large-halo" protogalaxies. "Small-halo" objects are 
an important fraction of the DM collapsed mass until redshift z ~ 5. Each curve has two thick 
portions: the thick section on the left shows "small-halo" objects, and the thick section on the 
right, clusters. In between, according to our definition, are "large-halo" objects. 
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produce H2 (Ferrara 1998), or destroy H2. In this paper, we focus on radiative feedback processes. 
From the first works on this subject (Haiman, Rees, & Loeb 1997; Haiman, Abel, & Rees 2000; 
Ciardi et al. 2000; Machacek, Bryan, & Abel 2001), the negative feedback from the background 
in the H2 Lyman- Werner bands was thought to be the main process that regulates the formation 
of "small-halo" objects. The background in the Lyman- Werner bands dissociates H2 through the 
two-step Solomon process, suppressing or delaying the formation of "small-halo" objects. But there 
are also positive feedback processes that can enhance the formation rate of H2 (Mac Low & Shull 
1986; Shapiro & Kang 1987; Haiman, Rees, & Loeb 1996; Ferrara 1998; Oh 2001; Ricotti, Gnedin, 
& Shull 2001). It is not trivial to determine which feedback prevails by means of semi-analytic 
models. Are "small-halo" objects able to form and survive the negative feedback effect of the disso- 
ciating background, or is their number drastically suppressed? Are "small-halo" objects numerous 
enough to have some effect on the subsequent evolution of the IGM? In this work we try to answer 
these questions, using self-consistent 3D cosmological simulations with radiative transfer. 

We present our results in two papers. This first paper focuses on the implementation of the 
code and its convergence. In Ricotti, Gnedin, & Shull (2002) (Paper II) we show and discuss the 
results. 

The simulations we present here are the first 3D cosmological simulations that include a self- 
consistent treatment of galaxy formation with radiative transfer by following the evolution of DM 
particles, star particles, gas, and photons. We also solve the line radiative transfer for the back- 
ground in the Lyman-Werner bands. A number of potentially relevant physical processes are 
included as well: secondary ionizations of H and He, detailed H2 chemistry and cooling processes, 
heating by Lya resonant scattering, H and He recombination lines, heavy element production and 
radiative cooling, and the spectral energy distribution (SED), g^, of the stellar sources consistent 
with the choice of the Lyman-continuum (Lyc) escape fraction, {fesc)i from the resolution element 
of the simulation. These ingredients, together with careful convergence studies, allow us to study 
the formation and feedback of "small-halo" objects by means of numerical simulations. The advan- 
tage of this approach, compared to semi-analytic models, is that we make no assumptions about the 
physical processes that are important in the simulation, apart from sub-grid physics. We assume 
that, when the gas sinks below the resolution limit of the simulation, it forms stars according to 
the Schmidt law (see eq. [20]). 

The paper is organized in the following manner. In § 2 we describe the cosmological code and 
the detailed physical processes simulated. In § 3 we study the convergence of the simulations, and 
in § 4 we summarize the features of the code and its limitations. We also provide a few examples 
of the results to demonstrate the power of the code to for addressing astrophysical problems. 
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2. The Code 

The simulations presented in this paper were performed with the "Softened Lagrangian Hy- 
drodynamics" (SLH-P^M) code described in detail in Gnedin (1995, 1996), Gnedin & Bertschinger 
(1996); Gnedin & Ostriker (1997), and Gnedin & Abel (2001). The code solves the system of time- 
dependent equations of motion of four main components: dark matter particles (P'^M algorithm); 
gas particles (quasi-Lagrangian deformable mesh using the SLH algorithm); "star-particles" formed 
using the Schmidt law in resolution elements that sink below the numerical resolution of the code; 
and photons solved self-consistently with the radiative transfer equation in the OTVET approxima- 
tion (Gnedin &; Abel 2001). We solve the line radiative transfer of the background radiation in the 
H2 Lyman- Werner bands with spectral resolution Av/v = 9.26 x 10""^ (i.e., 20,000 logarithmic bins 
in the energy range 11.2—13.6 cV of the H2 Lyman- Werner bands). Detailed atomic and molecular 
physics is included: secondary ionizations of H and He, H2 chemistry and cooling processes, heating 
by Lya resonant scattering, H and He recombination lines, heavy element production and radiative 
cooling, and galaxy spectrum consistent with the value of (fesc)- 

For clarity, we first summarize some relevant features of the code already implemented and 
published. The new components of the code are described in more detail as separate subsections. 

2.0.1. Radiative transfer: the OTVET approximation 

The OTVET approximation is based on solving equations for the first two moments of the 
photon distribution function. Since the moment equations have the conservative form, the photon 
number and flux conservation are achieved automatically in this approximation. However, the hi- 
erarchy of moment equations is not closed at any finite level, i.e. the flux (first moment) equation 
contains the second moment (fiux tensor), which can be reduced to the Variable Eddington tensor 
(the latter has a unit trace). In the OTVET approximation, the Eddington tensor is calculated 
in the optically thin approximation - hence the abbreviation OTVET: Optically Thin Variable 
Eddington Tensor approximation. Since the two moment equations can be solved very fast numer- 
ically, most of the computational time is spent in calculating the Eddington tensor, which can be 
done using the same algorithm used for computing gravity, in our case the Adaptive P^M. 

Methodically, this approach is equivalent to separating the algorithm for calculating absorption 
and emission of photons from the algorithm for calculating the direction of propagation of a photon. 
The first part can be done rapidly and accurately by using moments of the distribution function. 
The second part, in a direct ray-tracing approach, will be unacccptably slow. The OTVET ap- 
proximation allows us to speed up this calculation by an enormous factor without violating the 
accuracy of the first algorithm. 

Another feature of our implementation of the OTVET approximation is the use of "effective" 
frequencies. Accurate computation of photoionization and photoheating rates on a logarithmically 
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spaced mesh in the frequency space requires about 20 points per each e-folding, or about 300 

frequency bins altogether. Computing the radiative transfer for each of these bins is impractical at 
the moment. However, since wc arc ultimately only interested in the integrals over the frequencies, 
we can speed up the calculation of the reaction rates by a factor of about 50 by adopting the 
following ansatz for the radiation energy density 1,^ as a function of frequency: 



where index a runs over the list of species, which includes H I , He I , and He II (since we treat H2 in 
the optically thin approximation), Ii/^OT is the radiation energy density calculated in the optically 
thin regime, and the column densities, -/v|^\ are functions of position and time only. Thus, in 
order to compute Nj;^^ , we need to solve the radiative transfer equation in the optically thin regime 
(which is equivalent to just one frequency since the frequency dependence factors out) and at three 
other frequencies which we choose to be just above the Hi , Hel , and Hell ionization thresholds. 



In this version of the code we do not include the mechanical energy input from SN explosions, 
but we do account for the metal enrichment produced by SNe. The cooling function is calculated 
according to the gas metallicity. Since it is impractical to treat detailed ionization and thermal 
balance for some 130 ionization states of most common heavy elements (metals), we assume that 
whenever metal cooling is important, the different ionization states of heavy elements are in ioniza- 
tion equilibrium. In that case, we can adopt a single cooling function that describes the cooling rate 
due to excitation of atomic transitions in heavy elements. The metallicity of the stars produced 
in the simulation are recorded, but the emission spectrum is not calculated consistently with the 
metallicity of the stellar population. 



Wc adopt a ACDM cosmological model with the following parameters: CIq = 0.3, i}\ = 0.7, 
h = 0.7, and 0,h = 0.04. The initial conditions at z = 100, where we start our simulations, are 
computed using the COSMICS package (Bertschingcr 1995) assuming a spectrum of dark matter 
perturbations with power-law slope n = 1 and COBE normalization. We parameterize the unknown 
sub-grid physics of star formation at high redshift with three free parameters: the star formation 
efficiency, e*, the energy fraction in emitted ionizing photons per baryon converted into stars, euv, 
and the Lyman continuum (Lyc) escape fraction, {fesc)- 

In the rest of this section, we will describe in detail some new components included in the code 
that are particularly useful to address the problem of the formation of the first luminous objects. 




2.0.2. Heavy elements 



2.0.3. Cosmological model 
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2.1. Line Radiative Transfer in the Lyman- Werner Bands 

The evolution of the specific intensity Ji, [erg cm~^ s^^ Hz^^ sr~^ ] of ionizing or dissociating 
radiation in the expanding universe, with no scattering, is given by the following equation: 

^ + 4 i^J'^ SJ.) = -KJ. + S.. (1) 



dt dx^ ^ ' \ du 

Here, are the comoving coordinates, H is the Hubble constant, k^, is the absorption coefficient, 
S'jy is the source function, and x* = en'' /a, where n* is the unit vector in the direction of photon 
propagation and a = (1 + z)~^ is the scale factor. The volume-averaged mean specific intensity 
("background") is 

Mt) = {Mt,x,n))v, (2) 
where the averaging operator acting on a function f{x,n) of position and direction is defined as: 

{f{x,n))v = lim -L- [ d^x [ dnf{x,n). (3) 

The mean intensity Ju{t) satisfies the following equation: 

H(^^-3jJ]=-kJ, + S,, (4) 



dt \ d 

where, by definition. Si, = {Si,)v, and ki, = {k,jJy)v / Jv In general, ki, is not a space average of 
ki,, since it is weighted by the local value of the specific intensity J^. In the limit of the mean free 
path of radiation at frequency v being much larger than a characteristic scale one is interested in 
(in our case the size of a computational box), = {ki,)v- This limit is approached for the H2 
photodissociating radiation in the Lyman- Werner bands. 

If we rewrite equation (4) in terms of the dimensionless comoving photon number density 
riv = {aL})ox)^^T^Jv /hp where L},ox is a spatial scale (in this paper we take it to be the comoving 
size of the computational box), and hp is the Planck constant, using substitutions dt = {a^/Hojdr 
and ^ = ln(i/), equation (4) can be reduced to the following dimensionless equation: 

dnc 1 da dnc 

where = a^k^/HQ and = L^^^AirSi, / hpH^. Using a comoving logarithmic frequency variable, 

e = e + ln(a), (6) 

equation (5) can be reduced to 

drii 



dr 



« = -a^n^ + S^, (7) 
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which has the formal solution, 



n 



|-(r + At) = n|(r) exp 



+ 



/t+At r />T+Ar 

dt'S^{t')ex^p - J dr'a^ir') 



(8) 



We calculate equation (8) at each time step, At, of the simulation. The two integrals inside the 
square brackets on the right side of equation (8) can be solved analytically. We solve the third 
integral numerically. The analytical calculations are analogous to the ones presented by Ricotti 
et al. (2001) [see § 5 in that paper; in particular eqs. (28)-(31) and Fig. 11]. By definition, the line 
absorption coefficient is 



1 1 



Lyc 

i=LW i=Ly/3 



where cj) are the line profiles, fg^^ and /^^i ^'■^ oscillator strengths of the lines in the H2 
Lyman- Werner bands and Hi Lyman series lines, respectively, and f^yi^^ is the probability for the 
jth jjj^g^ calculated from Black & Dalgarno (1976), to decay to the ground vibrational level of the 
X^E+ ground electronic state of H2. A fraction fiyi=o of photons are re-emitted at the absorption 
frequency; therefore, these photons are not removed from the dissociating background. 

We assume that the absorption lines in Lyman- Werner bands have Gaussian line profiles, 
(j)G{y,yi) = (l/x/TrAi^j) exp[— (z/ — vi)^ / Avf], where Avi = 3 x 10~'^ ViXj^j^ is the Doppler width 
of the line i. The strong resonant lines in the hydrogen Lyman series typically have Lorentzian 
profiles (l)L{i^,Ui) = {ri/2n)/[{u - Vif + {Ti/2f], where = 7i/27r and 7^ = ^ is t^ie 

natural width of the z*'* H I line. Carrying out the integration, we have: 



-+A7 



dt'adt') 



vre 



fLW 

J osc,i 



Mv Yl ^(l-A."=o)^G(r*,Ar) + 



i=LW 
Lyc nH 

+{nm)v Yl -^^Lir*,AT) 



(9) 



for the lines in the H2 Lyman- Werner bands and the H I Lyman lines from Ly/3 upward. Note that 
Lya is not in the H2 dissociating band. The functions $g ^ind are given by 



erf 



Aui 



- erf 



arctan 



Aui 
— arctan 



F,: 



(10) 
(11) 



Here, v' = u\\+H{t)At\, v" = z/[l— if(r)(r— r*)], and erf(x) is the error function. The temperature 
dependence of the Gaussian line profiles can safely be neglected. As shown in Figure 11 of Ricotti 
et al. (2001), the integrated line profile $g is almost independent of the gas temperature. In the 
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limit of TiGM = 0, the integrated line profile is the Heaviside function, 




1 if {ui - Ai < ly < i^i - Aji) 

if (z/ < — Aj) and [u > Ui — Aji). 



(12) 



The same method is used to calculate the H2 dissociation cross section a^is from the two-step 
Solomon process: 



where fdis,i is the fraction of absorptions in the H2 line i that cascade to the dissociating continuum. 
In the top panel of Figure 2 we show the first and the second terms on the right-hand side of 
equation (8). In the bottom panel we show the dissociation cross section cr^is (dashed line) and the 
source term (solid line). In this figure we assume At = 1 Myr, z = 30, and XH2 = 10~^. The 
spectral resolution is Av/v = 9.26 x 10~^. 



In Figure 3 we compare the chemical rates from the three compilations used most widely in 
cosmology. The labels on the top of each panel indicate the reaction and the name of the rate 
coefficient according to the author of the compilation: K (Shapiro & Kang 1987), A (Abel et al. 
1997) and H (Galli & Palla 1998). The bottom-right panel compares the ro- vibrational H2 cooling 
functions from Tegmark et al. (1997), Lepp k Shull (1983), and GalU & Palla (1998). We have 
assumed xh,^ = 10^^; for comparison we also show the Lya cooling. The thin lines show additional 
cooling from the reactions k7, k9, kl6, and kl7 according to Shapiro & Kang (1987), assuming 
Xg = 0.5,Xjj+ = 10~^ and Xjj- = 10~^. These cooling processes are generally negligible. The 
cooling from reactions A8 and AlO, Aj^j = nH{3.53ksnH- + 1.38/cion^+) eV cm~^ s~^ (Abel et al. 
1997), can be important (horizontal dash-dotted line). The code can use each one of these rates, 
but in the following we will use the rates and H2 cooling function from Galli & Palla (1998). We also 
include the additional cooling/heating processes from the aforementioned reactions. We neglect D 
and Li reactions and the formation of HeH"*" molecules. They are negligible in the range of densities 
resolved by our simulations. 



Generally, in scattering processes of photons with particles, the recoil of the particles causes 
the radiation field to exchange energy with the matter. In the case of resonance lines, contrary to 
electron scattering, the cumulative effects of atom recoils remain small even for large line optical 
depths. This happens because a photon in the red wing of a line is strongly biased to scatter back 




(13) 



2.2. H2 Chemistry and Cooling/Heating Processes 



2.2.1. Heating by resonant scattering of Lya photons 
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Fig. 2. — The three functions (in dimensionless units, relative to one) used to solve line radiative 
transfer in the H2 Lyman- Werner bands and Hi Lyman series. Top panel: first (solid line) and 
second (dashed-line) terms on the right hand side of equation (8). Bottom panel: the source 
term (solid line) and H2 dissociation cross section adis (dashed line). The spectral resolution is 
Au/i' = 9.26 X 10~^ (i-e., 20,000 logarithmic bins in the energy range 11.2 — 13.6 eV). 
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Fig. 3. — Comparison of rate compilations found in literature: Shapiro & Kang (1987) dashed lines, 
Galli &: Palla (1998) dash-dotted lines, and Abel et al. (1997) solid lines. The bottom-right panel 
shows H2 cooling functions from Tegmark et al. (1997) (solid line), Lepp & Shull (1983) (dashed 
line), and Galli &: Palla (1998) (dot-dashed line). The thin lines in the same panel show cooling 
rates from the reactions, k7, k9, kl6, kl7 (Shapiro k Kang 1987) and A8, AlO (Abel et al. 1997). 
In the simulations shown in this paper, we adopt the rates and the H2 cooling function from Galli 
& Palla (1998). We also include cooling from the aforementioned reactions. 
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to the blue. Thus, the background intensity at the resonance hne frequency develops only a slight 

asymmetry. The average relative shift in a Lya photon energy, E, after having been scattered by a 
hydrogen atom at rest is {AE/E) = —{hpu/mHC^) ~ —10^*^, where niH is the mass of the hydrogen 
atom (Madau, Meiksin, &; Rees 1997) and hp is the Planck constant. If T <C hpVa/k = 1.2 x 10^ K, 
the energy is transferred from Lya photon to the gas at a rate, 



where. 




is the Lya scattering rate per H atom and ctj/^" the Lya absorption cross section. 

The background at the Lya frequency is produced mainly by redshifted non-ionizing UV 
photons (i.e., the photons emitted at the rest frame in the Lyman- Werner bands). All the other 
sources of Lya photons have a negligible effect on heating of neutral IGM. For instance, compared 
to X-ray heating discussed in the next section, the heating produced by the scattering of Lya 
photons emitted in the rest frame (not redshifted) from normal galaxies, is about 100 times less 
efficient (Madau et al. 1997). 



2. 2. 2. Secondary Ionization and Heating from X-rays 

Photoionization of H I , He l , and He II by X-rays and EUV photons produces energetic pho- 
toelectrons that can excite and ionize atoms before their energy is thermalized. This effect can be 
important (Venkatesan, Giroux, &; ShuU 2001) before reionization, when the gas is almost neutral 
and the spectrum of the background radiation is hard due to the large optical depth of the IGM 
to UV photons. 

We provide analytic fits to the Monte Carlo results of Shull & van Steenberg (1985). Collisional 
ionization and excitation of Hell by primary electrons arc neglected since, in a predominantly neu- 
tral medium, they are unimportant. The primary ionization rate for the species i = H I, He I, He II 
are, 

C = 47r/ -^eM-r.Xdv, (16) 

where Ti, is the continuum optical depth, and a'l is the photoionization cross section of the species 
i. Secondary ionizations enhance the photoionization rates as follows: 

C"' = C"^+ E C(^"'(^0,^e)) (17) 

i=Hi,Hci,Hcii 

^Hc^^He,^ E C(^"^^(^0,^e)), (18) 
i=Hi,Hei,He li 
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where ($"'(E^,Xe)) and {^^^'{Elxe)) express the average number of secondary ionizations per 
primary electron of energy Eq = hpV — P weighted by the function (Jj^/ZipZ/) exp(— r,/)(T* . Here 
/* = hpVi is the ionization potential for the species i. 

The photoionization heating rates for the species i = H I, He I are given by, 

r = 47r / eM-r)(TlEh(El,, Xe)dv. (19) 

In Appendix A we give analytic fits to the functions <I>^', ^^^^ and E^. We express these functions 
in a form suited to minimize the computational time for the integrals in equations (17)-(19). 



2.3. Sources spectrum and (/esc) 

In Ricotti et al. (2001) we showed that the SED of a galaxy and (fesc) are the main parameters 
that determine the relevance of positive feedback effects. In particular, we found that the Popula- 
tion HI (metal-free) and mini-quasar (with spectral index a = 1.8) SED produce similar positive 
feedback regions (PER). In this section, we explain how these ingredients are included in the code. 
We do not use the mini-quasar SED in our simulations. Apart from producing a stronger X-ray 
background, their effects should be analogous to those of the Population HI SED. 

Resolution elements that sink below the resolution limit of the simulation are allowed to form 
stars and therefore become sources of radiation according to the following equations: 

Star formation is implemented using the Schmidt law, where is the star formation efficiency, 

and pg are the stellar and gas density, respectively, and is the maximum of the dynamical time 
and cooling time. The parameter ejyy is the ratio of energy density of the ionizing radiation field 
to the gas rest-mass energy density converted into stars, [photons Hz~^] is the normalized SED, 



f 



g^dv = 1, 

'I/O 

and T,y is the internal opacity of the source. The escape fraction of ionizing photons, {fesc): is 
defined as: 

/•oo 

if esc) = / gue-^-diy, (22) 

J I/O 

where vq is the H i Lyc frequency. A value of (fesc) < 1 has the additional effect of producing harder 
source spectra because the gas optical depth is higher for UV and FUV photons than for X-ray 
photons. We include this effect modifying the SED using a frequency dependent optical depth 

T,. = -/Vhi(<7hi + aoo-Hei + aio-Heii)- (23) 
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Here, oq = -/Vnei/A^Hi and ai = Nuen/^Hi where Ni is the column density of the species/ion i. The 
unknown parameters, oq and ai, depend on the interstellar medium (ISM) properties of the galaxies, 
including clumping. Based on the results of simple ID radiative transfer simulations presented in 
Ricotti et al. (2001), we chose ao = 0.01 and ai = 10 when we use the Population III SED, and 
ao = 0.1 and ai = 10 when we use the Population II SED. These values are calculated placing the 
source of radiation in the center of a spherical galaxy with a gas density profile calculated with 
high-resolution numerical simulations (Abel, Bryan, & Norman 2000). We did not consider the 
effects of a clumpy ISM. Given the value of {fesc), -^Hi is calculated from equations (22)-(23) using 
a bisection algorithm. 

The spatial resolution of the code is determined by the dimensionless SLH parameter, S*. 
The code switches from Lagrangian to Eulerian when the deformation of the spatial grid exceeds 
a critical value that is approximatively 1/-B* the undeformed cell size. Therefore, the comoving 
spatial resolution is about Li,ox/{NboxB^), where Lfeoi is the box size and iV^^^ is the number of 
cells in the box. Usually we use = 10 for N^ox = 64 (262, 144 cells), = 16 for Ni,ox = 128 
(2 million cells) and B^ = 25 for Nbox = 256 (17 million cells). The gravity force is calculated 
by a P'^M algorithm with the Plummer softening parameter equal to 1/2B^. This ensures that 
the gravity force is accurate on our resolution scale, and thus 1/-B* is a real rather than nominal 
resolution. The nominal resolution is twice higher. 

We allow the spectral energy distribution, gi,, to be chosen from the following: 

i) Population II: metallicity Z = 0.04 Zq, evolutionary tracks evolved to t = 1 Gyr, continuous SF 
law, Salpeter IMF with star masses between 1 < M* < 100 M0 (Leitherer et al. 1995). Wolf-Rayet 
stars are responsible for the substantial EUV emission in this SED. 

ii) Population III : metallicity Z = 0, non-evolved, instantaneous SF law, Salpeter IMF with star 
masses between 1 < < 100 M0 (Tumlinson &; Shull 2000). Note the importance of Hell ionizing 

radiation in this SED. 

iii) Population III : metallicity Z = 0, evolutionary tracks evolved to i = 1 Gyr, continuous SF law, 
Salpeter IMF with star masses between 1 < M* < 100 M© (Brocato et al. 1999, 2000). This SED 
calculation differs substantially from (ii) since it does not include He II ionizing photons. 

In Figure 4 we show the analytical fits implemented in the code for these three SEDs. In our 
simulations we will use only the non-evolved Population III SED (ii) and the Population II SED 
(i). The evolved Population III SED (iii) is implemented in the code but never used. We show it 
just for comparison since calculations from different authors of the Population III SEDs are still in 
disagreement. 

The free parameter ejjv is the ratio of energy density of the ionizing radiation field to the gas 
rest-mass energy density converted into stars (/9*c^). Thus, ejjv = {hpi^ /mHC^)(^y-, where e^y 
is the number of ionizing photons per baryon converted into stars. The value of ejjv depends on 
the IMF and the metallicity of the stellar population. Assuming a Salpeter IMF with star masses 
between 1 Mq < M* < 100 M©, we have {euv/'^'n:) = 1.1 x 10-^ 75^/(13.6 eV) = 1.76 for the 
Population II SED and (e[/y/47r) = 2.5 x 10-^ 71^/(13.6 eV) = 2.47 for the (ii) Population III 



DRAFT: February 1, 2008 



15 



10 



■14 



10 



■15 



10 



16 _ 



N 



10 



18 



10 



■19 



10 



-20 



"1 I I I I I I I 



"1 I I I I I I I 




J I I I I I I I 



J I I I I I I I I 




10^ 



10< 



hpiy[eV] 



Fig. 4. — Analytical fits to the normalized Population II SED for a continuous SF law evolved to 
t = 1 Gyr (solid line), non-evolved Population III SED for an instantaneous SF law (short-dashed 
line) and Population III SED for a continuous SF law evolved to i = 1 Gyr (long-dashed line). See 
text for further explanation. 
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SED. Theoretical works (Uehara et al. 1996; Larson 1998; Nakamura & Umemura 1999, 2001) 

have shown that the IMF of metal-free stars could be dominated by massive stars or might have a 
bimodal IMF (lacking in intermediate mass stars). If the IMF is dominated by very massive stars 
(M ~ 300 M0) that eventually evolve into black holes (Bromm, Kudritzki, &: Loeb 2001), it is 
possible to have (e[/y/47r) > 2.5 x 10~^. 

The free parameter {fesc) is the fraction of ionizing photons escaping from the resolution 
element of the simulation. According to its definition, this parameter is, in general, resolution 
and time dependent. Theoretical studies (Ricotti & ShuU 2000; Wood Sz Loeb 2000) have shown 
that (fesc), defined as the fraction of ionizing photons escaping from a galaxy halo-*^, decreases at 
increasing redshift and increasing mass of the halos (assuming constant star formation efficiency). 
In particular, Ricotti & Shull (2000) provide analytic formulae for (fesc) {zvir,^, fg,MuM) as a 
function of the star formation efficiency (SFE), e, normalized to the Milky Way value, the redshift 
of virialization, z^ij., the DM halo mass, Mjjm, and the fraction of collapsed gas fg. The following 
formula, for example. 



is valid for "small-halo" objects {Mdm < W^Mq[{1 + Zvir)/10]'^-^), a constant SFE 4 < e < 400 
(i.e., the ionizing flux is proportional to the gas mass of the halo) and assuming a power law 
for the ionizing photon luminosity function of the OB associations with slope a = 2 (Dove & 
Shull 1994; Kennicutt, Edgar, & Hodge 1989), lower limit Si = 10^^ photons s~-^, and no upper 
limit. Equation (24) cannot be used in the code, because the definition of (fesc) depends on the 
spatial resolution of the simulation. Instead, (fesc) in equation (24) is a function of the halo mass. 
Nevertheless, even if (fesc) in the code is larger than (fesc) in equation (24), depending on the spatial 
resolution of the code, the functional dependence of (fesc) on the virialization redshift should be 
robust. In this study, we will consider {fesc) as a constant free parameter. Equation (24) is simply 
used to provide a first-order estimate of reasonable values for {fesc) at high redshift. We plan to 
relax this assumption in future work. 



In this section we investigate the numerical convergence of our simulations. Simulating the 
formation of the first galaxies is computationally challenging. It requires high mass resolution in 
order to resolve the internal structure of each low-mass object with a sufficient number of dark 
matter particles and grid cells. Moreover, since the first objects arise from the collapse of the 
most massive and rare density perturbations at redshift z ~ 30, we need a sufficiently large box to 



^The comoving radius at 200 times the background density of a DM halo with mass Mdm that virializes at redshift 
Zyir is raoo w (2.4 h- "^kpc) {Mo M/W^h-^MQy/^Kl + z)/{l + Zyir)]. 




(24) 



3. 



Resolution Studies 
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include at least a few 3 a perturbations of the initial density field. These requirements translate 
into the need for simulations with large mass dynamical range (i.e., a large number of cells and 
DM particles in the box). We have been able to run a few simulations with the maximum feasible 
resolution of 256^ cells, but the bulk of the study consists of runs with 128^, and 64^ cells. 

In Table 1 we list the runs used to test the convergence of the simulation. In all these test 
runs (but one: 256Llp3) the radiative transfer algorithm has been shut down. Therefore, the stars 
were able to form but did not produce any feedback effects. 



3.1. Mass resolution and Box Size 

We check the convergence of the simulations with two sets of runs. The first set that includes 
runs 64L1, 128L1, and 256L1 has a constant box size L^ox = 1 h"^ Mpc, and mass resolutions 
Mdm = 3.15 X 10^,3.94 x 10^, and 4.93 x 10^ h'^ Mq, respectively. The second set that includes 
runs 64L05, 128L1, and 256L2 has a constant mass resolution Mjjm = 3.94 x 10'^ Mq and 
box sizes L^ox = 0.5, 1, and 2 Mpc, respectively. If we assume that a physical output of the 
simulation, C, scales with N^ox as 

(a usual assumption in convergence analysis of numerical simulations) , we can determine the asymp- 
totic limit Coo = ^^^Niox^ooC if we have three simulations that differ only in the value of N^ox- 
Since the choice of equation (25) is somewhat arbitrary, the estimate of the convergence of the 
simulations is perhaps not rigorous and the results should be interpreted simply as a sign of con- 
vergence. We never used resolution studies to correct the results our the simulations for resolution 
or box size effects. 

The curves in Figure 5 show equation (25) as a function of Njjgx, where the observable C is the 
redshift z at which the star fraction is log f star = log{pstar / Pb) = const. In Figure 5 (left) the data 
points are for the first set of simulations; therefore, the asymptotic limits Coo (horizontal lines) show 
the rcdshifts at which log fstar = —6, —5, —4, —3, —2.2 for a simulation with infinite mass resolution 
and Lhox = 1- Figure 5 (right) is analogous to Figure 5 (left) but for the second set of simulations. 
Here, the asymptotic limits Coo show the redshifts at which log fstar = —6, —5, —4, —3.5 for a 
simulation with infinite box size and mass resolution Mdm = 3.94 x 10^ Mq. 

In Figure 6 we summarize the mass resolution and box-size convergence studies. The curves 
show the fraction of baryons in stars, fstar, for the two sets of simulations discussed in the previous 
paragraph. The triangles show fstar in the limit of a simulation with infinite mass resolution 
and Lf,ox = 1 Mpc. The squares show fstar for a simulation with infinite box size and mass 
resolution Mdm = 3.94 x 10^ Mq. Figure 6 shows that the 256L1 run is the closest to the 
convergence limit. 
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Table 1. List of simulations. 



RUN 


^box 


{h-^ Mpc) 


Mass Resolution 
{h-^ Mo) 




( ) if esc) 


e* 


64L05B8 


64 


0.5 


3.94 X 10^ 


8 





0.2 


64L05 


64 


0.5 


3.94 X 10"^ 


10 





0.2 


64L05B12 


64 


0.5 


3 94 X 10^ 


12 





0.2 


64L05B16 


64 


0.5 


3.94 X 10^ 


16 





0.2 


64L1 


64 


1.0 


3.15 X 10^ 


10 





0.2 


64L1B16 


64 


1.0 


3.15 X 10^ 


16 





0.2 


128L05 


128 


0.5 


4.93 X 10^ 


10 





0.2 


128L1 


128 


1.0 


3.94 X 10^ 


10 





0.2 


128L1B16 


128 


1.0 


3.94 X 10^ 


16 





0.2 


256L1 


256 


1.0 


4.93 X 10^ 


10 





0.2 


256L2 


256 


2.0 


3.94 X 10"^ 


10 





0.2 


256Llp3^ 


256 


1.0 


4.93 X 10^ 


25 


2.5 X 10-6 


0.1 



Note. — Parameter description. Numerical parameters: Nj^^^ is the number of 
grid cells, -L^ox is the box size in comoving h^^ Mpc, is a parameter that regu- 
lates the maximum deformation of the Lagrangian mesh: the spatial resolution is 
~ Lhox/i^boxB*). Physical parameters: is the normalized SED, e* is the star 
formation efficiency, euv is the ratio of energy density of the ionizing radiation 
field to the gas rest-mass energy density converted into stars (depends on the 
IMF), and {fesc) is the escape fraction of ionizing photons from the resolution 
element. 

'^Qi, is the Population III SED, modified assuming {fesc) = 0.1, ao = 
Nuei/^Hi = 0.01 and ai = Nueii/Nui = 10 where iVj is the column density 
of the species/ion i (see § 2.3). 
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Fig. 5. — Resolution studies. (Left) we use a set of three simulations (256L1, 128L1, and 64L1) 
with different mass resolution and L^ox = 1 Mpc, to extrapolate to the infinite resolution limit 
(thin horizontal lines). The curves, from top to bottom, show z + 1 as a function of N^ox [eq. (25)] 
at constant log /stars = —6, —5, —4, —3, —2.2, and crosses on each curve show the three data points. 
(Right) same as left but for a set of three simulations (256L2, 128L1, and 64L05) with different 
box sizes and the same mass resolution. The thin horizontal lines show the infinite box size limits. 
The lines from top to bottom are at constant log f stars = —6, —5, —4, —3.5. 
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Fig. 6. — Fraction of baryons in stars versus redshift for two sets of simulations with different 
resolution and box sizes. The three thick lines show 256L1, 128L1, and 64L1 simulations with 
constant box size {Li,ox = 1 Mpc), and varying mass resolutions {Mdm = 3.15 x 10^, 3.94 x 10^, 
and 4.93 x 10'^ h^^ M©, respectively). The triangles show fstar as a function of redshift in the limit 
of a simulation with infinite mass resolution and Lhox = 1 h^^ Mpc. The three thin lines show 
64L05, 128L1, and 256L2 simulations, with constant mass resolution {Mdm = 3.94 x lO'' h'^ Mq) 
and varying box sizes {Lhox = 0.5,1, and 2 Mpc, respectively). The squares show fstar as a 
function of redshift in the limit of infinite box size and mass resolution Mdm = 3.94 x 10^ /i"^ M©. 
Note that the simulation 128L1 appears in both sets; therefore there are only five lines in the plot. 
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The level of convergence of a simulation can be easily understood from the plot in Fig- 
ure 7. The cusp-shaped curves show the fraction of cooled gas in collapsed objects at 2; + 1 = 
35, 30, 25, 20, 15, 10, 4 (from bottom to top) calculated with the Press-Schechter formalism. The 
gas collapsed in the DM halo is assumed to be shock-heated to about its virial temperature. In 
order to form stars, the gas cooling time tcooi{z) must be shorter than the Hubble time^ 1/H(z). 
We roughly estimate the fraction of baryons able to cool and form stars by multiplying the viri- 
alized object fraction by the cooling efficiency ecooi = ''nax[l/tcooi{z)H{z),l]. Note that ecool < 1 
for masses typical of both "small-halo" objects and galaxy clusters. The comoving size of the box 
determines the DM mass of the most massive halo in the simulation, and the mass resolution de- 
termines the mass of the least massive object resolved. The oblique lines show the DM mass of the 
most massive halo in simulations of box sizes Ljjo^ = 0.5, 1, 2 Mpc (from top to bottom). The 
vertical lines show the mass of the least massive DM halo for which we resolve the SF in simulations 
with mass resolution Mdm = 4.93 x 10^,3.94 x 10"^, and 3.15 x 10^ /i"^ M© (from left to right). 
As a quick reference, the vertical lines correspond to 256^, 128^,64'^ cubes with L^o^ = 1 h^^ Mpc. 
These minimum masses arc about 100 times the mass resolution of the simulation. That means 
that we need to resolve each halo with about 100 cells. 

Given the box size and the mass resolution of the simulation, the point of intersection of the 
appropriate vertical and oblique lines shows the redshift of formation and the mass of the first 
object in the simulation. For example, the two thicker vertical and oblique lines delimit the masses 
resolved in the 256^ cells, L^o^ = 1 Mpc simulation. In this simulation, the first objects have 
masses M ~ 5 x 10^ Mq and form between 2: = 29 and z = 34. Moreover, this simulation can 
resolve the bulk of the collapsed baryons up to redshift z ^ b. 

In this study we show both integrated physical quantities such as SFR, background radiation 
intensity, chemical abundances, and also the properties of individual bound objects formed in the 
simulations. We identify bound objects with the DENMAX algorithm (Bcrtschinger & Gelb 1991). 
The results of DENMAX depend on the choice of the free parameter, G = L^io^/Lg-, which defines 
the smoothing length, Lq, of the particle density field. In Appendix B we show how to select the 
smoothing parameter when comparing simulations with different resolution/box sizes. 

In Figure 8 we compare the SFR at 2; = 15 as a function of halo mass for simulations with 
different resolutions and box sizes. The halos have been identified using DENMAX with smoothing 
parameter G = 400 for the 128 boxes and G = 1000 for the 256 box. It is evident that the 128L1 
run does not sufficiently resolve small mass {Mdm < 2 x 10^ M©) halos, and the 128L0.5 run 
does not have enough large mass objects. 

The minimum redshift to which we can evolve the simulations depends on the box size. For 2, 
1, and 0.5 Mpc boxes, we can evolve the simulation until redshifts z ^ 6,9, and 12 respectively. 



^For a ACDM cosmology with {n^, Ha, h) = (0.3, 0.7, 0.7) the age of the universe at a + 1 = 35, 30, 25, 20, 15, 10 
and 4, is tn = 82, 103, 136, 190, 292, 537 Myr and 2 Gyr, respectively 
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Fig. 7. — Fraction of collapsed gas as a function of DM halo mass at z + 1 = 35, 30, 25, 20, 15, 10, 4 
(from bottom to top) . We use the Press-Schechter formalism and a simple estimate for the fraction 
of collapsed gas based on the cooling time (see text). The thin portions of the curves correspond 
to virial temperatures 10^ ^ T^r ~ 10^ K. In these halos, the cooling is very efficient and all the 
gas is collapsed; the thick portions of the curves on the left are "small-halo" objects {%ir < 10"^ K) 
and on the right objects with T^r > 10^ K. The vertical lines correspond to the smallest halo mass 
fully resolved by simulations with mass resolution Mdm = 4.93 x 10^,3.94 x lO'^, and 3.15 x 10^ 
Mq (from left to right). As a quick reference they correspond to 256^,128^,64^ cubes with 
Lbox = 1 Mpc. The oblique lines show the largest halo mass that we can find in cubes with 
Lbox = 0.5, 1, 2 Mpc. For example, the two thicker vertical and oblique lines delimit the galaxies 
fully resolved in the 256L1 simulation. 
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Fig. 8. — Global SFR at z = 15 as a function of DM halo mass. The three histograms show the 
same simulation at different resolutions and box sizes. The solid line shows the simulation with 
256^ cells and Lhox = 1 Mpc, the dashed-hne 128^, Lhox = 1 Mpc, and the dotted hne 128^, 
Lbox = 0.5 h^^ Mpc. The halos have been identified using DENMAX with smoothing parameter 
G = 400 for the 128^ boxes and G = 1000 for the 256'^ box. The three horizontal lines above the 
histograms show the total (integrated over all masses) SFR for the same simulations. 
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Thereafter, scales of the order of the box size become nonhnear, and the box ceases to be a 
representative sample of the universe. 

3.2. Spatial resolution 

The spatial resolution of the grid is regulated by the softening parameter B^,, which sets the 
maximum deformation of tlic grid (the cell size can only become ~ 25* times smaller than the initial 
value) . The value of cannot be arbitrarily large because the grid has a "numerical tension" that 
prevents excessive deformation. 

In Figure 9 wc show the fraction of baryons in stars, fstar^ as a function of redshift for values 
of = 8, 10, 12, 16. It appears that the results are not sensitive to the value of B^, at least in the 
explored range. The only noticeable effect is a slight delay in the formation of the first objects as 
B^ increases. 

4. Discussion and Summary 

We have introduced a 3D cosmological code with radiative transfer suited to simulate the 

formation of the first objects self-consistently. We use the "Softened Lagrangian Hydrodynamics" 
(SLH-P^M) code (Gnedin 1995, 1996; Gnedin & Bertschinger 1996; Gnedin & Abel 2001) to solve 
the system of time-dependent equations of motion for dark matter particles (P^M algorithm), 
gas particles (quasi-Lagrangian deformable mesh using the SLH algorithm), and "stellar-particles" 
formed using the Schmidt law in resolution elements that sink below the numerical resolution of the 
code. The highest mass resolution that we achieve is Mdm = 4.93 x 10^ M©, and the highest 
spatial resolution is 156 pc comoving. We solve the continuum radiative transfer with the 
OTVET approximation (Gnedin & Abel 2001), and wc solve exactly the line radiative transfer in 
the H2 Lyman- Werner bands of the background radiation. Detailed atomic and molecular physics 
is included, although we have not included "supernova feedback" (work in progress). 

Including ionization from secondary electrons in the simulation does not produce significant 
changes of the IGM temperature or ionization of H and He. Secondary electrons are important when 
the ionizing spectrum is hard and UV photons do not dominate the ionization rate. X-rays (keV 
energies) in the background radiation are produced either by Wolf-Rayet stars in Population II 
or by the much hotter Population III stars. Nevertheless, the ionizing radiation background is 
dominated by the UV photons in the spectral features produced by the redshifted resonant He I 
and Hell Lya lines. In Figure 10 we show three examples of background (i.e., spatially averaged) 
radiation spectra: for a Population III SED and (fesc) = 1, for a Population III SED and (/esc) 
= 10~^, and for a Population II SED and (fesc) = 10~^. The spectrum is rather complex, but we 
can identify spectral fcatTircs associated with the redshifted resonant Hi , He I , and He 11 Lya lines, 
and the jumps at the Hi , Hel , and Hell ionization thresholds. In the spectral region between 11.2 
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Fig. 9. — Fraction of baryons in stars versus redshift for simulations with different values of the 
spatial resolution parameter, B^. The two thick lines have N^ox = 128, Li,ox = 1 Mpc {B^ = 
10,16). The thin solid and dotted lines have A^hox = 64, L^ox = 1 Mpc (i?* = 10,16). The 
other thin lines have N^ox = 64, Lbox = 0.5 h^^ Mpc {B^ = 8, 10, 12, 16). The spatial resolution is 
about Lbox/iNboxB*). 
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- 13.6 eV, responsible for the photodissociation of H2, we solve the line radiative transfer using a 
grid with 20, 000 frequency bins. An enlarged figure of this spectral region would show more clearly 
the fine absorption features associated with H2 and resonant Hi line opacity in the IGM. For the 
rest of the spectrum radiative transfer is solved for both the background and the spatially varying 
radiation dividing the spectrum in 4 frequency intervals [(i) optically thin between 1 - 11.2 eV, (ii) 
Hi ionizing radiation, (iii) He I ionizing radiation and (iv) Hell ionizing radiation] resolved with 
300 logarithmic bins altogether. 

In Figure 11 we show an example of 3D rendering of Hll regions at 2; = 27.6 from one of 
our simulations (iV(,oj; = 128, L^o^ = 1 Mpc) with radiative transfer. The small Hll regions 
produced by "small-halo" objects (the first "large-halo" objects form much later, cooling by Hi 
Lya, Sit z ^ 18) are clustered along the dense filaments of the cosmic large scale structure. In 
Figure 12 we show the global SFR as a function of redshift in the 256Llp3 run (thick solid line) 
compared to the same simulation without radiative feedback effects (thin solid line) and to a case 
where we do not allow the formation of any "small-halo" objects (thin dashed line). It is evident 
that the bursting mode of SF appears in the simulation with radiative feedback. The formation of 
"small-halo" objects is not dramatically suppressed by radiative feedback effects. In Paper II we 
study the nature of the feedback mechanism and the dependence of the SFR on the free parameters 
((/esc)) e*, ejjv, and g^) related to the unknown sub-grid physics (stellar IMF, star formation 
efficiency, and ISM properties at high-redshift) . 

Two important processes are not yet included in the simulations: H2 self-shielding and feedback 
from SN explosions. In some of our simulations with radiative transfer, we find that H2 in the 
filaments has a column density Nh^ > 10^^ cm~^, sufficient to shield the outgoing radiation emitted 
by embedded sources and the ingoing radiation emitted by external sources. We could crudely 
include the effect of H2 self-shielding for the outgoing radiation emitted by each source embedded 
in the filaments, modifying the source spectrum in the Lyman- Werner bands according to the ID 
radiative transfer calculations shown in Figure 5 in Ricotti et al. (2001). A complete treatment of 
line radiative transfer would necessarily require some approximation, since it is impossible to achieve 
the frequency resolution needed for the exact calculation. We have decided not to implement any 
of these two possibilities, based on the first results of the simulations. In Paper II we will show 
that the SF is not suppressed by the dissociating radiation, even if we neglect H2 self-shielding. 
The inclusion of self-shielding will not affect the global SFR in a significant way, since the primary 
feedback process that regulates the star formation is not H2 photodissociation. 

Supernova explosions could be important sources of feedback, both of heat and ionization. 
They could also produce a self-regulating global star formation and contribute to the transport of 
metals in the low density IGM. Unfortunately, we do not yet understand the correct modeling of 
their dynamical and thermal effects on the ISM/IGM. An approximate treatment of SN mechanical 
and thermal energy input is currently implemented in the code. We did not use this feature of the 
code, because we prefer to study the effects of radiative feedback and mechanical feedback from 
SN explosions separately. Moreover, we are not confident that the treatment of SN explosions in 
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Fig. 10. — Spectrum of the background radiation, J^. We show three representative examples: 
(i) using a Population III SED with {fesc) = 1 (solid line); (ii) using a Population III SED with 
if esc) = 10^^ (dashed line); (iii) using a Population II SED with (fesc) = 10~^ (dotted line). The 
dotted vertical line shows the He i ionization frequency, and the dashed vertical line shows the He ii 
ionization frequency. 




Fig. 11. — Example of 3D rendering of Hil regions at z = 27.6. The colors code the logarithm 
of H II abundance, "small-halo" objects produce small H ii regions clustered along the large scale 
structure filaments. 
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Fig. 12. — Global SFR versus redshift. The thick solid line shows the 256Llp3 run (Nfjox = 256, 
Lbox = 1 Mpc, Population III SED, {fesc) = 0.1, = 0.1, and (e,7V'/47r) = 2.5 x 10"^) 
that includes continuum and line radiative transfer and all the physics discussed in this paper. 
As a comparison, we show the same simulation without radiative transfer (thin solid line) and a 
simulation {Nhox = 64, L^ox = 1 Mpc) that forms only "large-halo" objects since we did not 
include H2 cooling and radiative transfer (thin dashed line). 
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the code is appropriate, given the complexity of ISM in galaxies. The investigation of the effects 
of mechanical feedback from SNe explosions will be the subject of our future work. 



This work was supported by the Theoretical Astrophysics program at the University of Col- 
orado (NASA grant NAG 5- 7262). The simulations presented in this paper were performed using 
SGI/CRAY Origin 2000 array at the National Center for Supercomputing Applications (NCSA). 



A. Secondary Ionization and Heating from X-rays 

For photoelectron energies Eq 3> 100 eV we use the fits given by Shull & van Steenberg (1985). 
For lower energy photoelectrons we have made fits to the Monte Carlo results shown in Figure 3 of 
the aforementioned paper. The functional form of the fits allow us to integrate in photon frequency, 
factoring out the dependency on Xg. We pre-calculate the integrals and store the results in three- 
dimensional tables for each frequency bin in order to avoid computationally intensive calculations. 

The fraction of secondary ionizations per primary electron of energy Eq is ^^{EQ,Xe), where 
Eq = hpV — P with P ground-state ionization potential of the species i and 

^^{Elxe) = [yi{j)fi{Eh) - y2{j)f2{Ei,)], (Al) 

where 

fo ifS^<28eV fo if£;^<28eV 

^ \ 1 if > 28 eV ^ \ (28 e\/Elf-^ if E^ > 28 eV. 

Here, since we neglect ionization and excitation of Hell from energetic primary electrons, i = 
Hi, He I, He II and j = Hi, He I. Since a fraction of the energy of primary photoelectrons goes to 
secondary ionizations, the heating rates from H i and He i photoionizations are less efficient by the 
complementary factor. 



Eh{El,Xe) = E'oil - yi{heat)fi{heat)+y2{heat)f2{heat)], (A3) 



where 



(0 ifSA<lleV (0 if^^<lleV 

h = i f2 = < (A4) 

[l if > 11 eV [(11 eV/Elf-'^ if E^ > 11 eV. 

The high energy limits of equations (A1)-(A3) are the functions yi that are given by Shull & van 
Steenberg (1985), and y2 has the form, 

yi = Cil-xlr (A5) 
y2 = Cx'^il-xlr, (A6) 

with < Xe < 1. The coefficients for equation (A6) are listed in Table 2 and the fits are shown in 
Figure 13. 
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Table 2. Fitting Coefficients. 



Parameter 


C 


a 


b 


c 


yi(Hi) 


0.3908 




0.4092 


1.7592 


yi(Hei) 


0.0554 




0.4614 


1.6660 


yi{heat) 


1.0 




0.2663 


1.3163 


2/2(Hl) 


0.6941 


0.2 


0.38 


2.0 


2/2 (He i) 


0.0984 


0.2 


0.38 


2.0 


y2{heat) 


3.9811 


0.4 


0.34 


2.0 




Fig. 13. — (Left) Fraction of energy deposited as Hi ionization per primary electron of energy Eq, 
as a function of the electron fraction Xg- The data points and curves are for Eq = 28, 50, 100, 200 
eV from bottom to top. (Right) Fraction of energy deposited as heat per primary electron of energy 
Eq, as a function of x^- The data points and curves are for Eq = 2, 11, 28, 50, 200 eV from top to 
bottom. 
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B. Group-finding Algorithm 

The DENMAX algorithm (Bertschinger &: Gelb 1991) identifies halos as maxima of the smoothed 
density field. The algorithm assigns each particle in the simulation to a group (halo) by moving 
each particle along the gradient of the density field until it reaches a local maximum. Unbound 
particles are then removed from the group. 

The results of DENMAX depend on the degree of smoothing used to define the density field. 
A finer resolution in the density field will split large groups into smaller subunits and vice versa. 
This arbitrariness in results is a common problem of any group-finding algorithm, and it is not 
only a numerical problem but often a real physical ambiguity. Especially at high redshift, since the 
merger rate is high, it is difficult to identify or define a single galaxy halo. 

Here we compare the mass function obtained using DENMAX with different smoothing scales 
Lq = Lbox/G with the Press-Schechter formalism. The main aim is to choose the parameter G 
consistently in order to compare the results from simulations with an arbitrary number of cells in 
the box N^^^. The mass of a DM particle is Mp = 8.26 x 10^°[(L6o^/l/i-iMpc)3/Ar3^^(J2^/0.3)] 

Mq and the comoving space resolution Lres = -^6ox/(-B*iVjoa;)- 

1/3 

If we impose that each group has a minimum of Np DM particles, we have Lq > Np' Lres, 
which implies, 

G<^^<B^N,ox, (Bl) 

where is the resolution parameter. The inequality on the right-hand side is obtained by imposing 
Np > 1. Another constraint on G can be obtained using a smoothing length comparable to the 
halo virial radius Lq ~ Rvir ~ 3.4 kpc(M/lO^M0)^/^ of the smallest object that we want to find: 

Note that the mass of the halo is M = NpMj-es- Combining equations (Bl) and (B2), we get 

^^J^ « 6.75 < B,. (B3) 

'box 



In Figure 14 we show the mass function compared to the Press-Schechter formalism with a 
top-hat window function. Better agreement with the Press-Schechter formula is obtained using 
G (4 — 5)7Vbox [i.e., using Np^ {2 — 5) in equation (B3)]. 
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Fig. 14. — Comparison of the mass functions (histograms) varying the smoothing parameter G and 
the Press-Schechter formaUsm with a top-hat filter (thick sohd hne). (Top-left) 64L05 and 64L1 
at z = 11 with G = 320,400. (Top-right) 64L05 at z = 19 with G = 80,160,320. (Bottom-left) 
128L05 and 128L1 at z = 12 with G = 320, 400. (Bottom-right) 128L1 with G = 400 and 256L1 
with G = 1000 at z = 15.13. 
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